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ABSTRACT 

With time-series spectroscopic observations taken with the Near Infrared Spectrom- 
eter (NIRSPEC) at Keck II, we investigated the atmosphere of the close orbiting 
transiting extrasolar giant planet, HD 189733b. In particular, we intended to mea- 
sure the dense absorption line forest around 2.3 micron, which is produced by carbon 
monoxide (CO). CO is expected to be present in the planetary atmosphere, although 
no detection of this molecule has been claimed yet. To identify the best suited data 
analysis method, we created artificial spectra of planetary atmospheres and analyzed 
them by three approaches found in the literature, the deconvolution method, data 
modeling via ^-minimization, and cross-correlation. As a result, we found that cross- 
correlation and x 2_ data modeling show systematically a higher sensitivity than the 
deconvolution method. We analyzed the NIRSPEC data with cross-correlation and 
detect CO absorption in the day-side spectrum of HD 189733b at the known plane- 
tary radial velocity semi-amplitude with 3.4<j confidence. 
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INTRODUCTION 



The study of exoplanets is one of the most vibrant and excit- 
ing fields in modern astronomy. In the past 18 years, more 
than 860 exoplanets have been discoverecQ. This has led 
to an increasing interest in the physical characterization of 
these new objects. As result of those characterization ef- 
forts, several chemicals have thus far been detected in the 
atmospheres of a few planets, as shown below. 

Most of the planets discovered so far are located too 
close to their host stars to appear separated. Consequently, 
the light coming from the planet and star are observed simul- 
taneously. When attempting to measure the light only from 
the planet, one has to find a strategy to remove the domi- 
nating star light. A prosperous strategy to learn more about 
the atmospheres of remote planets is through investigation 
of those planets that transit their parent stars. The atmo- 
spheric properties of transiting exoplanets can be measured 
in two ways: (1) during the passage of the planet in front 
of the star (transit), when part of the stellar light crosses 
the planetary atmosphere and the signature of chemicals 
in that atmosphere gets imprinted in the light we measure 
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from the star, and (2) during the passage of the planet be- 
hind the star (eclipse), when the star temporarily blocks 
the planet's emission and we can determine its tempera- 
ture, albedo, chemical comp osition, etc., from the diffe rence 
spectrum. Vi a transits, Nal (Charbonneau et al.ll2002l ). HI, 
OI and CII (Vidal-Madiar et al.ll200l l2004h .water vapor 
(|Tinetti et alfeOOTT ) have been detected from space in the at- 
mospheres of the h ot Jupiters HP 2 09358 b and HD 189733b. 
From the ground, Red field et alJ (|2008l) measured Nal in 
the a tmosphere of HD 189733b, while ISnellen et all (|2008l . 
120101 ) confirmed the Charbonneau et al. (2002) detection of 
Nal in HD 209458b and also detected CO via the analy- 
sis of the transmission spectrum using high-resolution spec- 
troscopy around 2.3 /im. The latter result allowed these au- 
thors to directly measu re the RV of a tran siting hot Jupiter 
for the very first time. iBean et al.l (|2010l ) investigated the 
atmosphere of the Super Earth GJ 1214b via transmission 
spectroscopy and found no evidence of atmospheric features, 
indicating a hydrogen atmos phere with high cloud s, or a wa- 
ter do minated atmosph ere dde Mooii et alj|2012h . Most re- 
cently, [Sing^et^al] (|201 lh and lColon et al.f i 2012h announced 
the first detections of potassium in two planets, XO-2b and 
HD 80606b. Furthermore, eclipses have already provided 
temperature measurements for over twenty planets, both 
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from space and more recently also from the ground (e.g. 
ISing fc L6pez-Moralesll2009l ; lLopez-Morales et al.ll2010T ). 

In this work, we focus on a further strategy to measure 
atmospheric features in planets that is based on intermedi- 
ate and high spectral resolution (i.e. R — A/AA > 20, 000; A 
denotes the wavelength) spectroscopy with very large tele- 
scopes. Key to this method is to observe a large number 
of spectral features coming from the planet and to observe 
them at different orbital phases so that the traveling faint 
planetary signal can be disentangled from the dominating 
stellar one. Considering this, the main advantage of this 
method is that it is not restricted to transiting planets (c.f. 
iBroei et al.ll2012l : iRodler et al]l2012h . 

In the past, this method was used in the optical with the 
goal to detect starlight reflected from hot Jupiters (i.e. mas- 
sive planets that are a few stellar radii away from their host 
stars) and to measure their exact masses. Although all those 
campaigns resulted in a non-detection of reflected starlight, 
stringent upper limits to the planet-to-star flux ratio and to 
the geometric albedo of these planets could be established. 
To date, the tightest 3<r-upper limits on the geometric albe- 
dos of the hot Jupiters r Boo b, HD 75289b and v And b 
are 0.39 (Leigh et al. 2003a), 0.46 (Rodler et al. 2008), and 
0.42 (Collier Cameron et al. 2002), respectively. These re- 
sults consequently provided important constraints on mod- 
els of the planetary atmospheres such as those by Marley 
et al. (1999) and Sudarsky et al. (2000, 2003). As a result, 
models that predicted a high reflectivity for the planetary 
atmosphere could be ruled out for the studied planets. 

Towards near-infrared (NIR; 1 — 2.5 /im) wavelengths, 
the planet-to-star flux ratios drastically i ncrease due to the 
strong thermal emissi o n of hot Jupiters. Wiedemann et al.l 
(l200lK : iDeming et all (l2005h; IBarnes et al.1 (|2007al lbl, |200S| . 
l20ld ) and ICubillos et al.l (|201 If ) observed hot Jupiters by 
means of high-resolution spectroscopy at near- infrared wave- 
lengths, but were not able to detect any molecules in their 
atmospheres. 

Very recently, our group (Rodler et al. 2012) as well 
as Brogi et al (2012) reported the first successful detec- 
tion of CO via high-resolution spectroscopic observations 
around 2.3 /im of the non-transiting hot Jupiter r Boo b. 
Both groups were able to measure the orbital motion of this 
planet, which allowed them to determine the previously un- 
known orbital inclination of the planet and to finally solve 
for the exact planetary mass. 

This paper is dedicated to the investigation of differ- 
ent data analysis approaches which were used by different 
research groups, with the goal to measure the planetary sig- 
nal via high-resolution spectroscopy (Sections 2 and 3). In 
the second part of the paper, we present our studies of the 
search for carbon-monoxide in the atmosphere of the hot 
Jupiter HD 189733b (Section 4) and a brief summary of our 
results and conclusive remarks (Section 5). 



2 METHODS 
2.1 Overview 

The idea of this approach is to observe spectral features 
in the planetary atmosphere via intermediate- and high- 
resolution spectroscopy. At visual wavelengths, the main 



contribution to the flux coming fr om a planet is the st arlight 
reflected from the companion (Seag er et al. 1998). This 
means that the planet reflects the stellar spectrum, which is 
shifted in wavelength with respect to the star due to the or- 
bital motion of the planet, and which is furthermore scaled 
down in intensity by a factor of several times 10 4 due to 
the albedo of the planet, the illuminated fraction of the 
planetary disk visible, and the size of planet. Towards in- 
frared wavelengths, the planet-to-star flux ratio drastically 
increases to the order of 10 -3 for hot Jupiters. At NIR wave- 
lengths, the planetary flux is almost entirely produced by 
thermal emission from the planet. To measure and identify 
atmospheric features in the p lanetary atmospheres, th eoret- 
ical models are required (e.g. Sh arp Sz Burrowsll2007l '). 

Key to this method is to observe a large number of 
spectral features in the planetary atmosphere to significantly 
overcome the planet-to-star flux ratio. Furthermore, it is im- 
portant to take a time-series of spectra and therefore to ob- 
serve the planetary features at different orbital phases of the 
planet, allowing to distinguish between the rather fixed stel- 
lar spectrum and the planetary one, which is periodically 
traveling in wavelength. 

Data analysis involves the subtraction of the domi- 
nating stellar spectrum as well as of the telluric spec- 
tru m of the Earth's atmo s phere at NIR wavelengths 
(see ICharbonneau et al.lll999l: ICollier Cameron et al1l2002l: 



Deming et al.l 120051 ; IBarnes et al.l l2007al . and IRodler et all 
20081. for detai le d des criptions). We note in passing that 
lLangford et al.l (|201 lh reports a different approach by 
searching for the planet and stellar spectra in Fourier space. 

By adopting one of the methods described in Sec- 
tions 12.21 - 12.41 the spectral features of the planet in each 
residual spectrum are then co-added to form a mean line 
profile of this spectrum. The following methods were used 
by different research groups: the deconvol u tion method , 



which was used b y ICo llier Camero n et al 



iLeigh et al.1 (|2003l l: IBarnes et al] l|2007al lbl. [200. 



2002); 



2010), 



and two straight-f orward data modeling app r oaches adopt- 
ing Y 2 -sta tistics (|Charbonneau et all 1 19991: Rodler et al] 
2008, 2010) as well as cross-correlat ion ((Snellen et al.ll2010l : 



Brogi et alj|2012l : IRodler et aT1l2012h . For each residual spec- 
trum, all three methods return - in the ideal case - a vector 
containing the mean line profile of the atmospheric features 
of the planet, which is Doppler shifted due to the instanta- 
neous RV of the planet at the time of the observations (cf. 
Fig.®. 

The next step involves the alignment of all mean line 
profiles in the time series with the alignment being a function 
of the RV semi-amplitude of the planet K p . To this end, we 
transfer each of the planetary line profiles from the velocity 
grid relative to the star (v) to a grid based on K p : 



sin 2-K(f> 



(1) 



where cf> is the orbital phase of the planet (0 = occurs at 
mid-transit for transiting planets), which is a priori known 
from the RV solution of the star. 

The total of the aligned mean line profiles of the time se- 
ries are then added up to form a single, overall mean line pro- 
file of the spectral features of the planet. We finally search 
for the global minimum or the maximum peak, respectively, 
for x 2 -statistics or cross-correlation / deconvolution method. 
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Once such a candidate signal at a specific K p has been found, 
the confidence level of the measurement is determined as de- 
scribed in Section 1231 

One of the most important aspects of this technique is 
that it allows the determination of the orbital inclination i 
of a non-transiting planet via K p and 



K p = -Kp.max sin i 



and 



/2ttG m, 

V -Porb 



1/3 



(2) 



(3) 



where G is the gravitational constant, m* the stellar mass, 
and P or b the orbital period of the planet. -Kp, ma x is the max- 
imum possible RV semi-amplitude that the planet can have, 
and which occurs for an orbital inclination i = 90°. The 
knowledge of the orbital inclination i allows us to finally 
solve for the true planetary mass m p = m Pjm i n / sini, where 
Tip.min denotes the minimum mass derived from the RV so- 
lution of the planet. We note that Eq. [3] is only valid for 
m* > m p . 



2.2 Deconvolution method 

The basic idea of this method is to deconvolve each ob- 
served star-free (and telluric-line-free) spectrum with a 
high-resolution reference spectrum that describes the atmo- 
spheric features in the planet, thereby s umming up all these 
features into one m ean line profile (e.g. iDonati et al.lfl997l ; 
iBarnes et aflll998t) . 

In general, the observed planetary spectrum s(x) at a 
certain position x on the detector can be approximated as 
a convolution of the template spectrum f(x') with the ap- 
parent mean planetary line profile p. The template spectrum 
f(x') can be an artificial reference spectrum of very high res- 
olution (R > 100, 000) or have the form of a list containing 
the positions and depths of the planetary absorption lines. 

The observed planetary spectrum is given by 



s(x) — / f(x') p(x — x) dx, 



(4) 



with j p(x) dx = 1. The discrete version of Eq. [4] is 



si= a 

k=i — n 



(5) 



where index i is the pixel number of the observed object 
spectrum, and k is the index of the num erical grid used for 
the intrinsic spectrum (|Endl et al.| [2000'l. n denotes a cut-off 
parameter of p with Pi-k = for \i — k\ > n, which allows 
us to shorten the infinite vector containing the mean line 
profile. 

Following the algorithm bv lEndl et all (|2000h . the grids 
of the reference spectrum / and the mean profile p can be 
oversampled with respect to the grid on which the observed 
object spectrum is recorded. The oversampled version of 
Eq. [5] follows as 



9(i+l)-l 

E 

j = qi 



( yi /*' Pj 



(6) 



where q — k/i is the oversampling factor, and j and k' the 
indices of the oversampled grids of the reference spectrum 
/ and the output vector p. We combine the terms and rear- 
range Eq. [6] as follows: 



- q g(i+l)-l- 



h' —q 



j = qi — k' 



fj) Pk> 



(7) 



where j is the index of the oversampled grid of the reference 
spectrum /, and k' is the index of the oversampled grid of 
the output vector p. Equation [7] is nothing but a matrix 
equation of the type 

-f = F -ft (8) 
with 

9-1 

FjV = f( q *i-k>+k)- (9) 

fe=0 

Let m be the dimension of the vector s, which repre- 
sents the observed object spectrum, and nq be the dimension 
of the oversampled vector ~p* containing the mean profile. 
Then the dimensions of the matrix F are m x nq. However, 
Eq.[H]is incomplete since each data point Sj of the observed 
data ~i exhibits its error Asj. The complete version of that 
equation follows as 

^ + aI = f^. (io) 

Now we are ready for the calculation of the output vector p 
containing the planetary signal by using a deconvolution al- 
gorithm; with this step, the signal-to-noise ratio (S/N) of the 
planetary signal can be boosted by a factor y/7 in the ideal 
case, where the I spectral features have the same weight. 

This problem constitutes an inversion problem which is 
ill-conditioned (due to As), but over-determined (the size 
of the object spectrum is much larger than the size of the 
mean line profile). There are several algorithms to mathe- 
matically solve for this problem. We tested different decon- 
volution approaches and found that least-squares deconvo- 
lution preserves best the planetary signal. Thus, we modify 
Eq. [10] to form a least squares problem: 



E 



(sj - F ifc gfc) 
(As,) 2 



(-f-F -tff E ("s>-F -f) -> min,(ll) 



where As is the error vector of "s^, E = Diag[Asj~ 2 , As^ 2 ] 
and m is the dimension of the vector of the observed object 
spectrum s. We find the least squares solution for the out- 
put vector containing the mean line profile of the planet ~p* 
by solving the matrix equation obtained by multiplying both 
sides of Equation [S] with F T E: 



F T E ~£ = F T E F -f 

^^ = (f t ef) 



(12) 
(13) 

The matrix resulting from the multiplication F T E F is 
square, symmetric and positive definite. Therefore, the in- 
verse matrix can be calculated by Cholesky decomposition 
(Press et al. 1992). 



2.3 Data modeling and x 2 -statistics 

For each of the residual spectra (i.e. stellar-free and telluric- 
free spectra), we create a planetary model M being a 
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theoretical spectrum fo r the planetary atmosphere / (e.g. 
ISharp fc Burrows! 120071 ) . which is degraded in spectral res- 
olution to the resolution of the observations. This version 
of the planetary model M is then Doppler-shifted as a 
function of K p sin 2ir<f) (remember that K p is the RV semi- 
amplitude of the planet, which is unknown for most of the 
non-transiting planets; cj> is the orbital phase of the planet 
which is a priori known from the RV solution of the star) 
with respect to the star and interpolated on the pixel grid 
of the observations. The spectrum is furthermore scaled in 
intensity as a function of a planet-to-star flux ratio (e\) at 
wavelength A as well as to the phase-function (i.e. the illumi- 
nation geometry of the observed planetary disk at different 
orbital phases - see e.g. iRodler et al]|2008r ). For pixel k, 

M k = e p f k {\ k (1 + A'psin27r0 c" 1 )}, (14) 

where c denotes the speed of light. 

Varying the free parameters K p and e, we finally search 
for the best-fit model M to all the residual spectra s by way 
of ^-minimization in appropriate search ranges, where the 
reduced \ 2 ls 

1 {M k - s k ) 2 . 



2 

Xv = 



N - 



(15) 



N is the total number of data points, and m is the number 
of fitted parameters. Note that As denotes the errors of s. 

2.4 Data modeling and cross-correlation 

For each residual spectrum s, we Doppler-shift the atmo- 
spheric model spectrum / of the planet by velocity v (see 
Eq. 1) in a given search range and interpolate it onto the 
pixel grid of s. We then calculate the normalized correla- 
tion degree C(v) for each r e sidua l spectrum following the 
formalism bv lCubillos etafl (|201ll ): 



C(v) 



W 



E {(/*(«)-/(»)) (**-»)} 
{e(/^)-/») 2 }{e(^-*) 2 } 



(16) 



where index k denotes the pixel number, / and s are the 
mean values of the atmospheric model spectrum and residual 
spectrum, respectively. The parameter w denotes the weight 
for the specific residual spectrum, while W is the sum of the 
weights of all residual spectra which are used in the cross- 
correlation. 



2.5 Determination of the confidence level 

Once a candidate signal of the planetary spectrum has been 
located, its confidence level needs to be determined. Note 
that this candidate feature is a peak in the case where 
the deconvolution method or the cross-correlation analysis 
has been employed. If the data was analyzed adopting x 2 - 
statistics, the candidate signal is reverted and produces a 
dip (Fig.©. 

One way to determine the confidence level of the candi- 
date signal is to employ a bootstrap randomization method 
(e.g. iKiirster et al.l I1997T ): random values of the orbital 
phases are assigned to the observed spectra, thereby cre- 
ating N different data sets (N = 100,000 in our analyses). 



Any signal present in the original data is then scrambled 
in these artificial data sets. For all these randomized data 
sets, we re-run the data analysis in given parameter search 
ranges and locate the candidate feature with its maximum 
value or its minimum xS.min value, depending on the em- 
ployed method. 

The confidence level of the candidate features is esti- 
mated to be w 1 - FAP = 1 - m/N, where FAP is the false 
alarm probability, m is the number of the best fit models 
having a xt value smaller or equal than the xt.min value 
found in the original, unscrambled data sets. We note that 
for the other two data analysis methods, we count the num- 
ber of the best fit models yielding an equal or higher peak 
than the peak value of the candidate signal. 



2.6 Determination of the error range of the result 

In the case of a significant detection of the planetary atmo- 
sphere (^ 3a, which corresponds to ^ 0.9973 confidence), 
we determine the la-error of the free parame ters K p and e 
via bootstrap resampling l|Barrow et al.|[l982h . To this end, 
we create randomly resampled data sets of the size of the 
original data set. In the resampling process, each set of spec- 
tra is randomly drawn from the general pool of spectra. For 
each spectrum a copy is kept for the artificial data set, but 
the original spectrum is returned to the pool with the effect 
that it can be drawn again. This way some of the original 
spectra will not appear in an artificial data set while oth- 
ers may appear more than once. Then, the parameters K p 
and e of the best fit model are stored. We create a total of 
10,000 artificial data sets to get a more precise estimate of 
the distribution of K p (and e). 

To determine the la-errors of the model parameters K p 
and e for the original data set we examine the distributions 
that were obtained for these parameters for the artificial 
data sets. We create 68.3% confidence intervals around the 
original K p and e values by taking the ranges adjacent to 
either side of the original values that each contain 34.15% 
of the artificial values. 



3 SIMULATIONS 
3.1 Data 

We created artificial data sets with the goal of testing 
the three data analysis methods and to find out which 
is most sensitive for the task. We adopted a PHOENIX 
model spectrum jHauschildt et alJI 19971 ; iHelling et aT]|2008l ; 
IWitte et alj|201lf) for a brown dwarf having a temperature 
of T = 1800 K and solar metallicity around 2.3 ^tm (Fig. [TJ 
lower spectrum). In this wavelength regime, the brown dwarf 
spectrum exhibits a dense forest of absorption lines almost 
entirely produced by the molecule carbon monoxide (CO). 

Each data set consisted of 79 spectra, which were 
Doppler-shifted as follows. We assigned an orbital phase 
value to each of the spectra, starting with (f> — 0.305 for 
spectrum 1 and subsequently increasing the value of (f> by 
0.005 for the following spectra. In total, the orbital phases 
ranged from <j> = 0.305 to 0.695 and therefore comprised the 
regions where a planet would appear at its brightest. Since 
hot Jupiters typically have RV semi-amplitudes of the order 
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Figure 1. Reference spectra adopted for the data analysis. The 
lower spectrum depicts a theoretical PHOENIX model for a 
brown dwarf having a temperature of T = 1800 K and solar 
metallicity. We used this spectrum to create the data sets as well 
as for data analysis. In the upper spectrum, the same spectrum is 
shown, but we removed 20% of the lines to study also the effect of 
missing lines in the data analysis (cf. Case 2 in the text). We note 
that this normalized spectrum has been shifted up for clarity. 



of K p — 100 km s _1 , we chose K p — 100 km s _1 and calcu- 
lated the value of the instantaneous RV shift v of the model 
spectrum by Eq. [T] 

We simulated the data for a cross-dispersed infrared 
spectrograph of high resolution. Thus, we degraded the spec- 
tral resolution of the spectrum to R = 60, 000 by a convo- 
lution of the model spectrum with a suitable Gaussian and 
interpolated the resulting spectrum onto a pixel grid rang- 
ing from 2.3 to 2.35 /im. This pixel grid consisted of a total 
of 4800 px, and one pixel corresponded to a velocity range of 
1.5 km s -1 . We then scaled the depths of absorption lines in 
the spectra due a chosen star-to-planet flux ratio e between 
10 -3 to 5 x 10~ 5 . In the final step, we added Poisson noise 
to the data in such a way that the S/N level was at 500 per 
spectral pixel, which is a typical value for high-resolution 
spectra in the infrared. 

We analyzed these artificial data sets by the three data 
analysis methods and investigated also the three following 
cases: In Case 1, we analyzed the data by adopting the 
PHOENIX model spectrum which we had used for creat- 
ing the data sets. For the analysis with the data modeling 
approaches, we further degraded the spectral resolution of 
the reference spectrum to 60,000. 

In Case 2, we explored the influence of missing lines in 
the reference spectrum, which was then adopted for the data 
analysis. The reference spectrum was a version of the brown 
dwarf spectrum, in which we had randomly removed 20% of 
the lines (Fig. [TJ upper spectrum). 

In the last case (Case 3), we studied how tiny continuum 
variations affected the result of the data analysis. To this 
end, we added a continuous sine wave to the modeled spectra 
having a period of 1000 pixel and a semi-amplitude of 10 -3 . 
Data analysis was carried out by adopting the PHOENIX 
model spectrum containing all absorption lines. 
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Figure 2. Output of the different data analysis methods for a 
single spectrum: deconvolution method (uppermost panel), data 
modeling with x 2 -statistics (middle) and cross-correlation (bot- 
tom). The data to be analyzed had been created in the way de- 
scribed in Section I3.ll by adopting a scaling value of e = 10~ 3 
and an RV semi-amplitude of K p = 100 km s _1 . 



3.2 Results 

For all three cases, we find that all three methods are al- 
most equally sensitive, with cross-correlation and \ -data 
modeling showing systematically a higher sensitivity than 
the deconvolution method (Figures [3] - . 

For Case 1 (CO spectra), we are able to retrieve the 
planetary signal at the correct value of K p at the 3a- 
confidence level down to planet-to-star flux ratios of e ~ 
1/8000 for all three methods. 

For Case 2 (fewer lines), we made use of the same data 
sets which had been created for Case 1, and analyzed them 
with the model spectrum for which we had randomly had 
removed 20% of the absorption lines. Fig. [4] illustrates that 
by employing cross-correlation, the data modeling approach 
using x 2 -statistics, and the deconvolution method, the plan- 
etary signal with a planet-to-star flux ratio of e « 1 /5500 is 
retrieved at the 3<r-confidence level. 

In Case 3 (variable continuum), we adopted the same 
data sets which we had created for Case 1, but multiplied 
the spectra with a sine wave to simulate variable continuum. 
Data analyses were carried out by adopting the correct spec- 
trum of the brown dwarf. As shown in Fig. [5] we are able to 
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Figure 3. Comparison between the data analysis methods for 
Case 1. We created different data sets for different planet-to-star 
flux ratios e (for clarity, we show the inverse of e) and analyzed 
them with the three data analysis methods. The solid, dashed 
and dotted graphs depict the false alarm probability (FAP) of 
the results obtained with the deconvolution method, the data 
modeling method using x 2 -statistics, and with cross-correlation, 
respectively. The four dotted horizontal lines mark the confidence 
levels in cr-units. As we compare these results, we find that all 
methods are almost equally sensitive, with cross-correlation and 
X 2 -data modeling showing systematically a higher sensitivity than 
the deconvolution method. 



retrieve the planetary signal at the correct position and at 
the 3(j-confidence level down to planet-to-star flux ratios of 
e ~ 1/7500 with all three data analysis methods. 

As we compare the individual methods and cases, we 
realize that the deconvolution method is systematically the 
least sensitive one of the three data analysis techniques. We 
attribute the ~ 3 — 7% lower sensitivity of the deconvolution 
approach to the additional processing layer, where a mean 
line profile is calculated at once according to a mathematical 
constraint (in our case: least-squares minimization). As a 
second point, for the deconvolution we adopt a reference 
spectrum consisting of delta functions (at the reference line 
positions) , which is sampled onto the same pixel grid as the 
observed spectrum. Due to that discrete sampling of the 
reference spectrum, the position of each reference line is at 
a full pixel, i.e. the center of the line is likely to be shifted 
by a fraction of a pixel. This, and the common case of line 
blending, where two lines might be treated as one, might 
produce distorted line profiles which then affect the overall 
mean line profile. 



4 REANALYSIS OF THE 

HD 189733B-NIRSPEC DATA SET 

4.1 The planet HD 189733b 

The parameters of HD 189733b and its parent star are 
provided in Table [T] This transiting planet is among the 
best studied exopl anets so far. Discovered in 2005 by 
iBouchv et al.l ([2005 ), it has quickly become the favorite tar- 
get for planet atmosphere studies, being located in one of the 
brightest known transiting system. Key res ults include the 
detection of a hot spot on the planet surface |Knutson et al.l 
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Figure 4. Same as in Fig. \3\ but for Case 2. We investigated 
how an incomplete spectral information in the reference spec- 
trum affects the data analysis. The data set was created with the 
correct brown dwarf model spectrum, while for the data analysis 
we adopted a version of that model spectrum, in which we had 
deleted 20% of the lines (cf. Fig. [TJ. In comparison to the other 
cases (Figures [3] and [5j , we see that the sensitivity of all three 
methods is strongly affected by an incomplete reference spectrum. 



< 



1e+00 



1e-01 



1e-02 



1e-03 



1e-04 




1a 



2a 



3a 



4a 



5000 6000 



7000 8000 9000 
star/planet ratio 



10000 11000 



Figure 5. Same as in Fig. [3] but for Case 3. We investigated 
how a variable continuum affects the data analysis. To this end, 
we added a sine wave to the data having a period of 1000 pixel 
and a semi-amplitude of 10 -3 . As result, we again find that the 
deconvolution method is significantly less sensitive than the other 
two data analysis methods. 



120071 . l2012h by studying the phase function of the planet 
in the infrared, indicating a temperature arou nd 1300 K, 
as w e ll as the discover y of high-altitude haze ()Pont et al.l 
120081 : ISing et all l201ll ). Several atoms and molecules in 
the p lanet atmosphere were found: water jGrillmair et al.l 
2008) , sodium seen in absorption at visual wavelengths 
(Rcd field et al.ll200St). as well as methane and carbon dioxide 
(e.g. ISwain et al.ll2009l: IPesert et alj|2009l : IWaldmann et alj 



l2012hTDesert et al.l (|2009| ) fou nd strong absorption around 
4.5 ja n, probably due to CO. iLecavelier Pes Etangs et aU 
(2010) measured strong evaporation of the planetary atmo- 
sphere due to the high irradiation from the host star. In 
addition to these discoveries, the spectrum of the dayside 
emission has been measured at NIR- and mid-infrared wave- 
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Table 1. 

planetary 
Bak6 
Bou5 
VV9 



Parameters of the star HD 189733 and 
companion. Abbreviations for the references 



its 



= iBakos et all i feOOrjl Bar7 
= iBouchv et all J2005h. Knu7 
Ivan Belle fc von Braunl 1120091 1 



Barnes et all ll2007a[) 
Knutson etall j2007li 



Parameter 


Value 


Error 


Ref. 


otar: 








Spectral type 


Kl-2 V 




Bou5 


K (mag) 


5.54 


0.02 


VV9 


m* (M Q ) 


0.82 


0.03 


Bou5 


R* (R©) 


0.76 


0.01 


Bou5 


T cff (K) 


5000 


50 


Bou5 


Planet: 








m p (Mj U p) 


1.15 


0.04 


Bou5 


Rp (Rjup) 


1.154 


0.033 


Bak6 


Teff (K) 


1300 


200 


Knu7 


a (AU) 


0.0313 


0.0004 


Bou5 


i (°) 


85.79 


0.24 


Bak6 


forb (d) 


2.218 5733 


0.000 0019 


Bak6 


t 0=o (BJD) 


245 3988.803362 


0.0023 


Bak6 


Kp (km s- 1 ) 


152.6 


2.0 


Bar7 



1.015 



1.010 



1.005 
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0.990 




2280 2300 2320 2340 
wavelength [nm] 



2360 



2380 



Figure 6. Upper spectrum: the CO model with spectral resolu- 
tion of 25,000. For clarity, the normalized spectrum is shifted up 
by 0.005 and scaled for a planet-to-star flux ratio of e = 10 — 2 , 
which is a factor 4.5 la rger than the actually measured value by 
IWaldmann et a.1.1 d2012t) . Lower spectrum: the two spectral orders 
used in the data analysis. The residual spectra are shown after 
the subtraction of the stellar- and telluric absorption lines as well 
as after the bad pixel removal. 



lengths (ICharbonneau et al.ll2008l ; lGrillmair et al.ll2008l ; and 
IWaldmann et al.ll2012l ~ 



4.2 NIRSPEC data and their reanalysis 

We r eanalyzed the data set published in iBarnes et all 
(2010), which had been taken with the goal to detect H2O 
and CO in the atmosphere of HD 189733b. Their data anal- 
ysis, however, resulted in a detection of low significance 
(98.8%) of these elements in the planetary atmospher e. Data 
were obtained with NIRSPEC (jlVlcLean et all 1 19981 ) at the 
Keck II Telescope, Hawaii, USA, on UT 2008 June 15 and 
June 22, when the dayside of the planet was almost en- 
tirely visible. A total of 373 spectra were recorded using a 
1024 x 1024 InSb Aladdin-3 array. The spectra were taken 
with a slit width of 0.432 arcsec, giving a spectral resolution 
of R w 25, 000. 

Using th e meth od outlined in lBarnes et all (|2007al ) and 
IBarnes et all |2008l ). we reduced the data and attempted to 
extract the planetary signature from time-series spectra by 
removing the dominant spectral contributions: namely the 
stellar spectrum and the telluric lines. Contrary to the data 
analysis of Barnes et al. (2010), we restricted the data analy- 
sis to the last two spectral orders comprising the wavelength 
region of A = 2.275 to 2.31 /zm and A = 2.347 to 2.383 ^m, 
respectively (Fig. [6]), where we expected the dense CO ab- 
sorption forest of the companion spectrum. We note that the 
latter spectral orde r was not used in the data analysis by 
IBarnes et all J2010T). In the wavelength regime of those two 
spectral orders. IWaldmann et all (| 20121 ) reported a planet- 
to-star flux ratio of e = 2.2 x 10~ 3 ~ 1/450 from secondary 
eclipse measurements of HD 189733b. 

The S/N of the residual spectra (i.e. after the removal 
of the stellar spectrum and telluric lines) was on average 370 
and 450 per spectral pixel in the first and second night, re- 
spectively. We rejected those frames taken when the planet 
was behind the star and not visible. In the end, we worked 



on 153 and 116 residual spectra, almost entirely covering 
the orbital phases <j> = 0.302 to 0.421 and = 0.515 to 
0.580, respectively for the first and second night, when the 
day-side of the planet was largely visible. We further iden- 
tified bad pixels and outliers and discarded them from the 
data analysis. We adopted cross-correlation (Section 2.4) to 
search for the CO spectrum of the planet in the residual 
spectra. As reference spectrum, we adopted a PHOENIX 
spectrum of a brown dwarf with a temperature of T=1500 K 
and with a spectral resolution of R = 25,000 (Fig. |6}. In 
the cross-correlation, we weighed the spectra according to 
their average SNR level per spectral pixel, and further ac- 
counted for both the systemic radial velocity of the star 
HD 189733 (v — —2.4 km s _1 ) and the barycentric velocity 
of the Earth. We then co-aligned and co-added the indi- 
vidual cross-correlation functions in the planet's rest-frame, 
thereby taking into account the orbital phase information of 
the planet at the barycentric Julian date of the observations. 



4.3 Results and discussion 

As shown in Figure [7] we find the strongest candidate fea- 
ture at K p = 154 km s _1 , which is located within the lcr- 
error range of the known RV semi-amplitude of HD 189733b 
being, K p = 152.6 ± 2 km s" 1 . 

Since the RV semi-a mplitude of the plane t K p had al- 
ready been estimated by IBarnes et all (|2007bl ) , we can re- 
strict the bootstrap randomization run to the search range 
K p = 152.6 ± 3a km s _1 (i.e. 146.6 to 158.6 km s" 1 ). As 
result of this bootstrap randomization run, we found that 
the candidate feature is at a confidence level of 99.54% and 
therefore represents a 2.9o"-detection of CO in the planetary 
atmosphere of HD 189733b. 

We furthermore carried out the data analysis on the two 
independent nights and found this candidate feature present 
in both nights (Fig. [7] lower panel). As a plausibility test, 
we varied the systemic RV of HD 189733 and determined 
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Figure 7. Upper panel: co-added cross-correlation functions in 
the rest-frame of the planet. The peak of the CCF occurs at an RV 
semi-amplitude of 154 km s _1 and appears close to the known RV 
semi-amplitude of HD 189733b being K p = 152.6 km s~ 1 (dotted 
vertical line). While bootstrap randomization runs showed that 
the peak of the CCF is at the confidence level of 99.54% (2.9a), a 
more straight-forward approach revealed that the candidate fea- 
ture is significant with 99.92% (3.4<r) confidence and therefore 
represents a detection of the CO absorption line forest in the plan- 
etary atmosphere spectrum of HD 189733b. Lower panel: same as 
above, but showing the results of the cross-correlation functions 
for the first night (solid line) and second night (dotted line). 



the confidence levels at the measured RV semi-amplitude of 
the planet, K p — 154 km s . The highest peak occurs at 
the genuine system velocity of HD 189733, which is v sys =- 
2.4 km s" 1 (Fig.©. 

In addition to that, we adopted a different strategy to 
determine the confidence level of the candidate feature. To 
this end, we calculated the cross-correlation values for all 
combinations of the two parameters K p and the systemic 
velocity « sys of HD 189733b, respectively, in the range of 
5 < K p < 200 km s" 1 and -100 ^ v sys sC 100 km s" 1 . The 
average value of the cross-correlation values for all those 
combinations is 0.00055, its scatter (rms) is 0.00368, while 
the peak value of the candidate feature shows 0.0129. We 
subtracted the average value from the peak value and di- 
vided the result by the scatter of the cross-correlation func- 
tions. This revealed that the candidate feature is indeed sig- 
nificant at the 99.92% (3.36<j)-confidence level. 



3 - 




-100 -75 -50 -25 25 50 75 100 
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Figure 8. As a test, we varied the systemic RV of the star 
HD 189733 and determined the confidence levels at the measured 
RV semi-amplitude of the planet, K p = 154 km s — 1 . The con- 
fidence levels had been determined by bootstrap randomization 
runs. A clear signal at the systemic velocity of HD 189733. being 
fsys=-2.4 km s — 1 (dotted line) is visible. 




relative RV to star [km/s] 

Figure 9. The individual cross-correlation functions (CCFs) of 
the 269 residual spectra of HD 189733b with the CO model spec- 
trum in the rest-frame of the star are shown. During the course of 
the observations, the orbital motion of the planet produces a RV- 
shift starting at about 150 km s _1 and ending at -75 km s — 1 , re- 
spectively for orbital phases of 0.30 and 0.58. However, no trace of 
the planet appears. The linear grey-scales indicate the strength of 
the cross-correlation signal (dark means absorption, bright means 
emission). 

When plotting the individual cross-correlation functions 
of the residual spectra with the CO model spectra, we do 
not see the planetary signature (Fig. |9} . Given the relatively 
large flux ratio between the dayside of HD 189733b and its 
host star, we should have been able to measure the trace 
of the planetary RV signal. Albeit the presence of system- 
atic noise that hamper the retrieval of the weak planetary 
signature, that may suggest low abundance of CO in the 
planetary atmosphere of HD 189733b. 

To demonstrate the power of this approach, we injected 
an artificial planetary signal into the residual spectra and 
analyzed these data sets with cross-correlation. We modeled 
the artificial planetary signal adopting CO spectra with solar 



© 2013 RAS, MNRAS OOO.ITHTTl 



Detection of CO absorption in the atmosphere of the hot Jupiter HD 189733b 9 




relative RV to star [km/s] 

Figure 10. We injected a spectrum showing CO with solar abun- 
dance, with a spectral resolution of R = 25, 000 and for planet- 
to-star flux ratio of 1/450 to the residual spectra of HD 189733b. 
The individual cross-correlation functions (CCFs) of the 269 spec- 
tra, to which we injected a planetary signal, with the CO model 
spectrum having solar abundance in the rest-frame of the star are 
shown. The linear grey-scales indicate the strength of the cross- 
correlation signal (dark means absorption, bright means emis- 
sion). It is shown that the trace of the RVs of the injected plane- 
tary signal can be recovered. 



abundance with a spectral resolution of R — 25, 000. The 
spectra were then scaled to an intensity ratio of 1/450 and 
shifted in RV according to the orbital motion of HD 189733b 
(i.e. for an RV semi-amplitude of K p = 153 km s _1 ). The 
results of the data analysis are shown in Figure 1101 The 
injected planetary signal could be recovered at the correct 
value of K p , and bootstrap randomization simulations for 
a search range of K p from 5 to 200 km s _1 revealed that 
the signal is significant at a confidence level of > 99.99. By 
adopting bootstrap resampling simulations, we furthermore 
determined the la-error ranges of the recovered planetary 
signal, being 2.5 km s _1 . Figs. [5] and [TU] also show some 
strong cross-correlation artifacts which are of the order of 
the injected planetary signal. We attribute these artifacts 
mainly to systematic errors coming from the removal of the 
stellar signal and the telluric spectrum. 

To estimate the flux ratio between the strong plane- 
tary lines Flines and the stellar continuum F 1 *, we again in- 
jected an artificial planetary signal into a scrambled set of 
the residual spectra and analyzed these data sets with cross- 
correlation. The injected planetary spectrum was scaled to 
a chosen intensity ratio Fii ne s/-F* and Doppler-shifted as 
described above. We recovered the injected the planetary 
signal and determined its confidence level by the aforemen- 
tioned straight-forward strategy. We found that for a scaling 
factor of Fknes/F* = 1.8 x 10 -4 , we attain 3.4er confidence. 
This value again indicates a low abundance of CO in the 
planetary atmosphere of HD 189733b, given the large flux 
ratio between the continua of the planet day-side spectrum 
and the stellar one, being e = 2.2 x 10 around 2.3 ^m. 



5 SUMMARY AND CONCLUSIONS 

We carried out studies to find out what data analysis ap- 
proach is best suited for the search for atoms and molecules 
in hot Jupiters via high-resolution spectroscopy. We first 
created artificial data sets consisting of spectra of planetary 
atmospheres, scaled them in intensity according to a chosen 
value of Flincs/F 1 *, and analyzed them by different data anal- 
yses approaches. As result, we found that the highest sensi- 
tivities to recover the weak planetary features are attained 
with cross-correlation and x 2_ data modeling, while the de- 
convolution method was less sensitive (~ 3 — 7%) than the 
two aforementioned methods. 

In light of these studies, we attempted to measure 
the dense CO absorption line forest around 2.3 micron 
in the day-side spectrum of the transiting hot Jupiter 
HD 189733b. By employing cross-correlation, we reanalyzed 
a time-series of spectra taken with the Near Infrared Spec- 
trometer (NIRSPEC) at Keck II during two nights and de- 
tect a candidate planet signal at an RV semi-amplitude 
K p — 154 km s , which is located within the lcr-error 
range of the known RV semi-amplitude of HD 189733b 
(K p = 152.6 ± 2 km s _1 ). While bootstrap randomization 
runs resulted in 2.9<r-confidence for this candidate feature, a 
more straight-forward test revealed that we detect the plan- 
etary signal with a S/N of 3.4. 

As a plausibility test, we independently carried out the 
data analysis for each of the two nights and found this can- 
didate feature clearly present in both nights. In addition, we 
varied the systemic RV of HD 189733 and found the highest 
peak at the genuine system velocity of HD 189733, which is 
D sys =-2.4 km s~^ 

In the past, iBarnes et ail l|201Ch analyzed the same data 
set in the wavelength region 2.21 - 2.36 /im and searched 
for water and CO absorption in the planetary atmosphere. 
These authors adopted a different data analysis strategy 
(different bad pixel correction, deconvolution method, dif- 
ferent spectral orders used) for their purposes and found 
in their data analysis a candidate feature of the planetary 
signal close to the RV semi-amplitude of the planet at the 
98.8% confidence level. 

We are consequently confident to claim a detection of 
CO absorption in the planetary atmosphere of HD 189733b. 
This work demonstrates the power of intermediate- 
resolution spectroscopy at infrared wavelengths to investi- 
gate the atmospheres of remote planets. The measured plan- 
etary CO signal is weak and may suggest a low abundance 
of CO in the planetary atmosphere of HD 189733b. In addi- 
tion, we observe CO in absorption, which indicates that the 
atmosphere of HD 189733b lacks a strong thermal inversion 
layer. 
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